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Magnetic confinement of the solar tachocline: 
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3D MHD global solar simulations coupling the turbulent convective zone and the radiative zone have been carried out. 
Essential features of the Sun such as differential rotation, meridional circulation and internal waves excitation are recov- 
ered. These realistic models are used to test the possibility of having the solar tachocline confined by a primordial inner 
magnetic field. We find that the initially confined magnetic fields we consider open into the convective envelope. Angular 
momentum is transported across the two zones by magnetic torques and stresses, establishing the so-called Ferarro's law 
of isorotation. In the parameter space studied, the confinement of the magnetic field by meridional circulation penetration 
fails, also implying the failure of the tachocline confinement by the magnetic field. Three-dimensional convective motions 
are proven responsible for the lack of magnetic field confinement. Those results are robust for the different magnetic field 
topologies considered, i.e. aligned or oblique dipole. 
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1 Introduction 



Since its discovery (iBrown et al.ll 19891) . the solar tachocline 
has puzzled the scientific community. In particular, a clear 
explanation for its extreme thinness (less tha n 5% of the so- 
lar radius Rq, see lCharbonneau et al is still lacking . 



The first theoretical work dedicated to the tachocline was 



carried out by lSpiegel and ZahnI (119921) . They considered an 
hydrodynamic tachocline and showed that such an interface 
layer will spread into the radiative interior because of ther- 
mal diffusion. As a result, they estimated that the differen- 
tial rotation should extend down to 0.3 Rq after 4.5 Gyears, 
which is in total contradiction with th e helioseismic inver- 
sions ( Schou et al.ll998l:[Thompson et al. 2003) . Consequently 
additional physical mechanisms have to be considered in or- 
der to properly explain the observed thinness of the tachocline. 

ISpiegel and ZahnI ( 1992 ) suggested that the stratifica- 
tion in the tachocline would imply an anisotropy in the tur- 
bulence, making it predominantly horizontal. It would then 
erode the latitudinal gradient of angular- velocity (lElliotJl997h . 
However .Gough and Mclntyre (1998) pointed out that such 
anisotropi c turbulence would on the con trary act as an anti- 
diffusion ( Dritschel and Mclntyrell2008l) . In fact, the ques- 
tion is still hardly settled when both radial and latitudinal 
shears of angular- velocity are taken into account jKiml2005l; 



tachocline would erase both the diffusive and anti-diffusive 
behaviors of 2D turbulence. 

iGough and Mclntyr3 ( 1998 ) then proposed that a (fos- 
sil) dipolar magnetic field confined in the solar interior could 
oppose the thermal spreading of the tachocline. This solu- 
tion offered also an explanation for the solid body rotation 
of the solar rad iation zone and such a fossil field was also 
invoked by iRudiger and Kitchatinov ([1997.) to confine the 
tachocline. The fossil field confinement scenario is in fact a 
double confinement problem. First, the imposed magnetic 
field erodes latitudinal gradients of angular-velocity, thus 
confining the tachocline. Second, the magnetic field has to 
remain confined in the radiation zone against its outward 
,ohmic diffusion. In order to confine the magnetic field, it 
was first argued that a meridional flow coming from the 
convection zone down to the radiative interior through the 
tachocline at the high latitudes could even tually prevent the 
field f rom diffusing o utward in this region ^Gough and Mclntyre 
19981) . At the equator. lwood and Mclntvrel(l201 ll) suggested 

that if the magnetic field we re to connect to the convection 

zone, the magnetic pumping ( Dorch and Nordlundl200ll: iTobias etal 



snears oi angular- velocity are taken into account (IKimlzuUM: 
Kim and LeDroyostl2007l: Leprovost and Kiml2006l: MieschI 



2003h . Further. iTobias et al.l (l2007l) showed that the intro- 
duction of (even weak) magnetic fields in the bulk of the 
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200 ll: Iziegler and Rudige3l2003k would 'confine' it below 
the region of intense shear. 

Despite the appeal of this scenario, none of the many 
numerical simulations that have been carried out in the past 

succeeded in recovering it comple tely. Th e first simulations 

of the radiation zone in 2D(GaraudMQ3) and in 3Dl Brun and ZahnI 
(l2006l) showed that if the confinement of the magnetic field 
fails, angular momentum is transported along the field lines 
into the radiation zone, making the radiative interior rotate 
differentially. The Ferraro's law of iso-rotation is then estab- 
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lished ( FerrarQ|[l937 ). in complete contradiction with solar 
observations. However these models did not allow for flows 
penetrating from the convection zone, su ch as plumes or 
meridi o nal circulation. In their simulatio n s. iGara ud and Garauc 
(l2008k : iRudiger and Kitchatinovl (l2007b : ISule et al.. (.20051) 
imposed a meridional circulation at the top of the radia- 
tion zone. They recovered partially the lGough and Mclntyre 
scenario, but one may object that those results are 
highly sensitive to the prescribed (profiles and amplitude) 
meridional circulation. Simulations coupling self-consistently 
the tw o zones of th e Sun were then carri ed out by RogersI 
(I2OIII) in 2D and bv lStrugarek et al.1 (l201lk (hereafter SBZl 1) 
in 3D in order to take into account the motions of the con- 
vection zone. Both studies were not able to validate the 
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ed flu 



bala 



magnetic confinement scenario. Finally. IWood et al.l (1201 1|) 
conducted an analytical study in a reduced cartesian model 
and proposed an i mproved theory for a magn etic confine- 
ment scenario a /f^ lOough and Mclntyrdd 1998k . 

We are interested here in the fact that SBZll showed 
that the convective motions were responsible for the lack of 
confinement of the buried magnetic field, and that in most 
of the works previously cited, the magnetic field enters the 
convection zone primarily at the equator. Since different 
magnetic topologies may lead to significantly different in- 
teractions with the convective motions, the universality of 
the results obtained with our axisymmetric dipolar field can 
be questioned. We thus conduct here numerical simulations 
based on the SBZll model to investigate the role of the 
magnetic field topology. In Sect. [2l we summarize the main 
ingredients of the model used by SBZll, and in Sect. [5] we 
study the impact of the magnetic topology on the confine- 
ment of the field. Conclusions and perspectives are reported 
in Sect.m 

2 Modeling the convection and the radiation 
zones 



FoUowinglBrun et al.l (120111) . SBZl 1, we use the ASH code 
dBrun et al.ll2004 to model 90% of the solar interior in 3D 
using the non-linear MHD equations under the anelastic ap- 
proximation. A LES (large eddy simulation) approach is 
used to parametrize turbulent diffusivity profiles; the de- 
tailed equations can be found in SBZl 1. 

Our initial hydrodynamical setup is exactly the same 
than the one used in SBZll. We recall here the principal 
features of the model. We use the solar rotation rate VLq = 
2.6 10"^5"\ the solar mass Mr, 



luminosity Lc^ = 3.8 10^^ 



= 2.0 10^^^ and the solar 
= 0.0 iu ^ erg . The background ther- 
modynamic q uantities were computed with the CES AM code 
(lMorellll997k . lence we use the solar stratification and ob- 
tain a realistic Brunt- Vaisala frequency N in the radiation 
zone. The radial gradient of entropy is negative above r^c^ = 
0.715 Rq, it defines the convection zone since the Rayleigh 
number (~ 10^) in this region is well above the critical 
Ra yleigh number for the onset o f the convective instabil- 
ity (iGilman and Glatzmaieij[T98ll) . We display in Fig. [T] the 
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Fig. 1 Radial energy flux balance normalized to the lumi- 
nosity of the Sun. The dashed vertical line at r = r^c^/^o 
denotes the base of the convection zone, the dotted verti- 
cal line at r = / Rq denotes the bounding radius of the 
initial magnetic field (see Eq. (O). 



radial energy flux balance normalized to L©. The energy 
is carried by the radiative flux in the radiation zone (below 
ri)cz^ the black dashed line), and the enthalpy flux transports 
the major part of the energy in the convection zone (above 

ncz)- 

The conv ective motion s exhibit banana-like shapes at 
the equator ( Br un and Too mre 2002), and more patchy pat- 
terns at higher latitudes. Angular momentum is transported 
mainly by convective Reynolds stresses in the bulk of the 
convection zone, establishing a differential rotation similar 
to the solar differential rotation (see for example the color 
contours in Fig.©. A tachocline develops, and the angular 
velocity shear extends down to r = 0.58 i?©. The merid- 
ional circulation which is self-consistently excited in the 
convection zone is roughly unicellular in each hemisphere 
(when time-averaged), with downflows near the poles and 
upflows at the equator. We stress here that those mean flows 
are self-consistently induced by the interaction between the 
convective motions, the rotation, and the baroclinicity in- 
duced by latitudinal gradients of entropy. Thus, they are not 
artificially imposed by parameters, by boundary conditions 
or prescribed profiles. The convective overshooting depth 
is dov = 0.04: Rq (see the enthalpy flux in Fig. [T]) and is 
known to scale as the square root of the filling factor of 
downflow plumes (IZahnlll99lh . Consequently, it is certainly 
overestimated in this model. The mean meridional circula- 
tion penetrates by due = 0.035 Rq below the base of the 
convection zone r^cz • Even if our model is not completely in 
the solar parameter regime due to our enhanced diffusivities 
(see SBZll), we observe a penetration of both the merid- 
ional circulation and the convection below r^c^ • We are thus 
confident that our model captures the key physical ingredi- 
ents of the upper tachocline. Finally, we observe also that 
the convective ove rshooting plumes excite gravity waves in 
the radiation zone dBrun etal.ll20TTl) . 
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Since we used enhanced viscosity thermal diffusivity 
n and magnetic diffusivity r], the time- scales involved in the 
simulation differ from what is occurring in the real Sun. Still 
we took great care to maintain the proper hierarchy between 
the time-scales, though they are somewhat closer to each 
other than what they are in reality. The ohmic time- scale 

= R%/r] is r^^ = 9.60 years in the convection zone, 
and = 192 years in the radiation zone. By comparison, 
the convective turnover time is Tconv = 28 days. We refer 
the reader to SBZl 1 for further details on the simulations. 



3 Influence of the magnetic topology 

3.1 On the axisymmetric dipole topology 

The expression of the magnetic field used in SBZl 1 is given 
by B = 5o {Br^r + Beee), with 



Br = — 

r- 

where ^(r, ( 
such that 



Bo 



^ sin 6 



(1) 



r sin 6A^n is constant on field lines and is 



^ = (r/Ry (r - Rbf sin^ forr<Rb, 
= for r>Rt, 



(2) 



where R^ = 0.57 Rq is the bounding radius of the confined 
field. One can immediately see that for such a choice of ^, 
the maximum radial gradient of Bg is located at the equator. 
Since the behavior of the magnetic field is mainly diffusive 
in the bulk of the radiation zone, it will evolve more rapidly 
where its gradients are maximum. As a result, it preferen- 
tially interacts with the tachocline (and eventually the con- 
vection zone) at the equator. 

One may construct modified 'dipoles' in order to move 
the location of maximum gradient with latitude. To do so, 
we consider a ^ function defined by 

* = [a sin^ (9 + (1 - a) sin^ cos^^ O] f{r) , (3) 

where a and p are related to the co-latitude 6max where 
the gradient of Bg is maximum, and /(r) controls the ra- 
dial shape of the dipole. We choose 6max = 10° and run 
the ASH code. We display in Fig. [2] the initial magnetic 
configuration, and the evolution of the magnetic field at a 
later instant. We observe that the dipole interacts with the 
convective motions primarily at latitude 80°, as expected. 
This confirms the key role played by the field topology in 
controlling the temporal evolution of the field in the stable 
radiative interior. 

Nevertheless, the end state in this case will be the same 
than in SBZl 1 . Even if a longer time is needed for the mag- 
netic field lines to interact with the convective motions at 
the equator, the modified dipole is still axisymmetric. As a 
result, its interaction with the convective zone will be simi- 
lar to the behavior described in SBZl 1 . Angular momentum 
will still be transported along the field lines, although slight 
differences may occur due to the little topological change in 




Fig. 2 Azimuthal averages of the angular velocity ft 
(color map) and the magnetic field lines (red lines), in the 
axisymmetric case. The ^} contours have also been tempo- 
rally averaged. 



the bulk of the radiative zone. This will result again in a dif- 
ferentially rotating radiative zone. In the end, the magnetic 
scenario for the solar tachocline will again certainly fail in 
this case. 

3.2 Non axisymmetric magnetic topology: the obhque 
dipole 

The oblique magnetic fields in stars have been a long stand- 

ing su b ject of research (iMestel and TakhaJl972l: l IMestel and WeissI 
1987h . For exam ple, A-type stars a re typically thought to be 
oblique rotators (IBrun et al.ll2005h . Given that a purely ax- 
isymmetric field establishes Ferraro's law, one may wonder 
whether an oblique dipole acts similarly. In addition, any 
confined dipolar magnetic field may explain the solid body 
rotation of the radiation zone. We consequently test the ro- 
bustness of the results of SBZll by considering a tilted 
dipole. 

We use the ASH formalism and decompose the mag- 
netic field into poloidal and toroidal components. 



V X V X (Cer) + V X (Aer) 



(4) 



The dipole aligned with the rotation axis is then simply writ- 
ten 

{ C = nr)Yo' ' 

where stands for the classical spherical harmonics. To 
tilt the dipole with an angle /3 with respect to the rotation 
axis, we simply write 



A = 

C = f{r) (cos/3ro' + 



(6) 
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Fig. 3 3D rendering of magnetic field lines. The volume 
rendering is the azimuthal velocity in the rotating frame in 
the south hemisphere, orange denotes positive v^p and blue 
negative . 

where /(r) = Bq {r/R^f (r - Ri,f and Ri, = 0.57 Rq is 
the bounding radius of the initial magnetic field. We choose 
/3 = 60° to tilt significantly the dipole. A 3D rendering of 
the magnetic field lines is displayed in Fig.O 

Although the magnetic topology is clearly different from 
SBZll, we observe that if (3 ^ 90°, then the initial mag- 
netic field is the sum of an axisymmetric dipole {i.e., similar 
to SBZl 1) and a non-axisymmetric component (see eq. ([S])). 
This is also made clear by comparing Figs. |4(a)| and |4(c)| 
which are taken at the initialization of the magnetic field. 
The azimuthally averaged field displayed in Fig. |4(a)| only 
retains the axisymmetric component of the tilted dipole, al- 
though the real magnetic configuration is oblique (Fig. |4(c)| ). 

Figs. |4(b)] is taken 70 convective turnover times later. 
We see that the axisymmetric component of the oblique 
dipole roughly evolves similarly to the pure axisymmetric 
dipole of SBZl 1 . We observe indeed that the initially buried 
magnetic field does not remain confined and connects the 
convection zone with the radiative interior. The non axisym- 
metric field likewise does not remain confined in Fig. |4(d)[ 
and its evolution is highly three dimensional. We stress that 
the field lines plotted in the radiation zone are a fair repre- 
sentation of the field topology, whereas in the tachocline and 
in the convection zone only the projection into the merid- 
ional plane of a 3D magnetic field with non zero azimuthal 
component is represented. 

The magnetic Reynolds number realized in the convec- 
tion zone is b elow the threshold for a dynamo to occur {Rm ^ 
100 vs 300 in lBrun et all (l2004l) ). As a result, the total mag- 
netic energy decreases as the simulation evolves. 

The polar slices in Fig. |4] seem to exhibit major changes 
in the Q profile at the poles, and we know from SBZll 
that the unconfined axisymmetric magnetic field will apply 




Fig. 4 Meridional cuts of the titled dipole case. The quan- 
tities in (a,b) are azimuthally averaged, while they are taken 
at the longitude (p ^ (i.e., the longitude where the pro- 
jection of the magnetic field on a polar slice is indeed a 
planar dipole, at least initially) in panels (c,d). The slight 
discrepancy between the dipole angle and 60° (or —30°) is 
due to the fact that we are not exactly Sit (p = 0. The color 
maps display the angular velocity ft, and the black lines the 
magnetic field lines. (a,c) are taken at the initialization of 
the magnetic field, and (b,d) approximately 70 convective 
turnover times later. 
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a torque leading to efficient transport of angular momentum. 
The evolution of the non-axisymmetric component may also 
have some effects on the angular momentum balance. Thus 
we examine the different terms of th e angular momen tum 
balance which are defined as follows (iBrun et al ] l2nn4h : 



North Pole 



dtipC) = -V-iF 



MC 



Hydro 



MHD 



(7) 

where the different terms correspond respectively to contri- 
butions from Meridional Circulation, Reynolds Stress, Vis- 
cous Diffusion, Maxwell Torque and Maxwell Stress. They 
are defined by 



r sin Op {{v'^v] 



sinO 



(8) 
(9) 

0) 



47r 
r sinO 



(5^)(BM)and 



(11) 
(12) 



where the subscript .m designates the meridional compo- 
nent of V and B. In the previous equations, we have de- 
composed the velocity and the magnetic field into an az- 
imuthally averaged part (.) and (/:?-dependent part (with a 
prime). The different contributions can further be separated 
between radial (J> along e^) and latitudinal (Te along eg) 
contributions. However here we only focus on the radial flux 
of angular momentum defined by 



Mr) 



sin 



dO 



(13) 



where {Oi , O2) maybe chosen to study a particular region of 
our simulation. 

We plot in Fig. [5] the time-averaged angular momen- 
tum balance near the north pole and near the equator about 
75 convective turnover times after the introduction of the 
magnetic field. The magnetic contributions are highlighted 
in red and are separated in axisymmetric (Maxwell torque 
MT) and non axisymmetric (Maxwell stress MS). At the 
beginning, those contributions are exactly zero. Since the 
field is initially buried deep in the radiation zone, we first 
notice that angular momentum is transported at the top of 
the radiative zone by (axisymmetric) Maxwell torque at the 
north pole (as already noticed in SBZll). This again leads 
to a differentially rotating radiative zone. What is remark- 
able here is that we also observe magnetic transport of an- 
gular momentum at the equator. The axisymmetric and non- 




Fig. 5 Time-averaged angular momentum balance (Ir) at 
the north pole {Oi = 60, 02 = 90, see eq. (fTSl) ) and at the 
equator (^1 = —15, O2 = 15). Xr is normalized to 10^^. 

4 Conclusions and Perspectives 



In a pr evious paper (SBZl 1), we tested the lGough and Mclntyre 

(Em 

magnetic confinement scenario of the solar tachocline 
with a 3D non-linear MHD model of 90% of the Sun which 
couples together the convective and radiative zones. We showed 
that a fossil axisymmetric magnetic field does not remain 
confined in the radiation zone, and thus cannot prevent the 
spread of the tachocline. 

We showed that the motions (both meridional circula- 
tion and convection) were not strong enough to confine the 
field and that upward motions at the base of the convection 
zone were on the contrary helping the field to expand inside 



the convection zone. Although our numerical experiment is 
axisymmetric (torque and stress) components contribute equalljfar from representing the real Sun, we are confident that our 
to the inward transport of angular momentum, thus acceler- parameter regime allows us to model the main ingredients 

of the tachocline confinement scenario. More precisely, we 
demonstrated in SBZl 1 that even if we consider high mag- 
netic diffusivity, our magnetic Reynolds number is higher 
than 1 in the overshooting/penetration region and in the con- 



ating the upper radiation zone at low latitude, while Maxwell 
torque slows down the high latitude. We can thus conclude 
that an oblique field is acting very similarly to a purely 
aligned dipole. 
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vection zone. The magnetic field evolution is thus domi- 
nated by ad vection and not by diffusion. The overshooting 
depth of the convective plumes is likely to be overestimated 
because our Peclet number is lower than the solar one, and 
the penetration o f the meridional circulation may be un- 



deres timated (see iGaraud and BrummelllbOOSl: IWood et al 



201 ih in our model. However, penetration and overshooting 



are clearly observed in our model. If more 'realistic' values 
of the parameters could have consequences on the polar dy- 
namics, we stress here that it may not have much influence 
at the equator, where the confinement of the magnetic field 
fails. 

The magnetic topology considered in SBZll was very 
simple (axisymmetric dipole), and the robustness of the re- 
sults is not obvious when considering the geometry of the 
field. We demonstrated in Sect. 13.11 that an axisymmetric 
dipole preferentially diffuses at the location where its max- 
imum radial gradient is, i.e., at the equator. This has some 
importance since SBZll showed that the lack of confine- 
ment speci fically at the equator was resp onsible for the fail- 
ure of the iGough and Mclntyra (Il998|) scenario. In order 
to test other topologies, we reported in Sect. 13.21 the study 
of an oblique dipole buried in the radiation zone. A tilted 
dipole can be decomposed into an axisymmetric and a non- 
axisymmetric part, and the axisymmetric component was 
shown to evolve similarly to the purely axisymmetric dipole 
of SBZll. As a result, the magnetic confinement also fails 
for the tilted dipole. 

The non-axisymmetric component of the titled dipole 
exhibits interesting three dimensional dynamics and trans- 
ports angular momentum at the equator. In order to get rid 
of the SBZl 1 dynamics and to isolate the effect of the non- 
axisymmetric dipole, simulations with a purely perpendic- 
ular dipole {i.e., with no axisymmetric component, /3 = 90 
in Eq. ©) will be reported in a future paper. 
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